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\ ABSTRACT 

■ We investigate the effect of dark energy on the density profiles of dark matter haloes 

' with a suite of cosmological N-body simulations and use our results to test analytic 

models. We consider constant equation of state models, and allow both w > — 1 and 
CNJ ■ w < — 1. Using five simulations with w ranging from —1.5 to —0.5, and with more 

| than ~ 1600 well-resolved haloes each, we show that the halo concentration model of 

, Bullock et al. (2001) accurately predicts the median concentrations of haloes over the 

range of w, halo masses, and rcdshifts that we are capable of probing. We find that 
the Bullock et al. (2001) model works best when halo masses and concentrations are 
. defined relative to an outer radius set by a cosmology-dependent virial overdensity. For 

CU a fixed power spectrum normalization and fixed-mass haloes, larger values of w lead to 

higher concentrations and higher halo central densities, both because collapse occurs 
£^ ■ earlier and because haloes have higher virial densities. While precise predictions of halo 

^ | densities are quite sensitive to various uncertainties, we make broad comparisons to 

Ci . galaxy rotation curve data. At fixed power spectrum normalization (fixed erg), w > —1 

quintessence models seem to exacerbate the central density problem relative to the 
standard w = —1 model. For example, models with w ~ —0.5 seem disfavored by the 
data, which can be matched only by allowing extremely low normalizations, as ^0.6. 
' Meanwhile w < — 1 models help to reduce the apparent discrepancy. We confirm that 

the Jenkins et al. (2001) halo mass function provides an excellent approximation to 
the abundance of haloes in our simulations and extend its region of validity to include 
models with w < — 1. 

Key words: cosmology: theory - dark matter - large-scale structure of universe - 
methods: N-body simulations 
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1 INTRODUCTION 

In the prevailing model of galaxy formation, galaxies assem- 
ble and evolve in the potential wells established by gravita- 
tionally bound haloes of cold and collisionless dark matter 
(CDM). Except for some possible difficulties on small scales, 
the CDM model is remarkably successful in explaining a 
large number of observations. However, this success requires 
an additional "dark energy" component, that drives an ac- 
celerated cosmic expansion, to be added to the universal 
energy budget. While the presence of dark energy is firmly 
established observationally, measuring its equation of state 
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as well as developing a theoretical understanding of the na- 
ture of the dark energy are two of the biggest outstanding 
problems in cosmology today. Dark energy not only affects 
the large-scale evolution of the Universe, but also the col- 
lapse histories and density structures of dark matter haloes. 
Understanding the precise nature of these effects is impor- 
tant for studies that aim to quantify the nature of dark en- 
ergy using strong (e.g., Sarbu, Rusin, & Ma 2001; Huterer 
& Ma 2004; Kuhlen, Keeton, & Madau 2004) and weak 
(e.g. Hu & Jain 2003; Bartelmann et al. 2002) gravitational 
lensing. Changing the dark energy model should similarly 
change expectations for galaxy rotation curves, and could 
affect one of the main small-scale problems facing CDM - 
the central density problem (e.g., Zentner & Bullock 2002; 
McGaugh, Barker, & de Blok 2003, and references therein). 
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In the present paper, we use a suite of N-body simulations to 
study how halo density profiles change as a function of dark 
energy equation of state, discuss our results in the context of 
analytic models, and discuss the observational implications 
of dark energy on galaxy scales. 

The existence of some form of dark energy is supported 
by a preponderance of data. Taken together, observations of 
the magnitude-redshift relation of type la supernovae (SNIa; 
Perlmutter et al. 1999; Riess et al. 2001; Knop et al. 2003; 
Barris et al. 2003), the power spectrum of cosmic microwave 
background (CMB) anisotropy (Spergel et al. 2003; Tegmark 
et al. 2003a), the power spectrum of galaxy clustering (Do- 
delson et al. 2003; Tegmark et al. 2003b), and the luminosity 
function and baryon fraction of clusters (Allen et al. 2003) 
provide nearly unimpeachable evidence for the existence of 
dark energy. The most common supposition is that the dark 
energy takes the form of a cosmological constant or vac- 
uum energy. In this case, the energy density p, and pressure 
p, are related through p — —p. An attractive alternative 
candidate for the dark energy is the potential energy of a 
slowly- varying scalar field (f>, or "quintessence" (e.g., Ratra 
& Peebles 1988; Caldwell, Dave, & Steinhardt 1998). 

A convenient parametrization of the dark energy is 
through an equation of state w = p^/p^ relating its en- 
ergy density and pressure. In general, the equation of state 
parameter w is time-varying, but it is useful to model 
quintessence with a constant equation of state parame- 
ter because current observational data sets have limited 
power to distinguish between a time-varying and constant 
equation of state (e.g., Kujat et al. 2002). Useful limits 
on the equation of state for the dark energy, assuming 
that it remains constant in time, come from SNIa studies, 
-1.67 < w < -0.62 (2cr; Knop et al. 2003), and can be 
refined by combining SNIa data with CMB anisotropy and 
galaxy clustering statistics yielding —1.33 < w < —0.79 at 
2a (Tegmark et al. 2003b). For our simulations, we adopt 
an empirical view and study five models with constant w, 
that span a comparably large range of parameter space: 
w = -1.5, -1.25, -1.0, -0.75, -0.5. 

Our theoretical understanding of halo profiles has im- 
proved recently largely through numerical simulations, per- 
formed in the context of CDM plus cosmological con- 
stant (ACDM) or standard CDM (SCDM, i.e. fi M = 1, 
fiq = 0) cosmologies (Navarro, Frenk, & White, 1995, 1996, 
1997, hereafter NFW; Kravtsov, Klypin, & Khokhlov 1997; 
Ghigna et al. 1998; Jing 2000; Bullock et al. 2001; Eke, 
Navarro, & Steinmetz 2001; Wechsler et al. 2002, hereafter 
W02; Zhao et al. 2003; Hayashi et al. 2003; Navarro et al. 
2003; for a review, see Primack 2003). It is generally un- 
derstood that the final density profiles of haloes are linked 
closely to their formation histories. Halo central densities 
are set during an early, rapid-accretion phase, and tend to 
be proportional to the density of haloes in the Universe at 
the time of this rapid collapse (W02; Zhao et al. 2003; Tasit- 
siomi et al. 2003). Larger values of w lead to earlier collapse 
times and also to more rapid collapse of overdensities, thus 
we expect haloes with higher relative central densities and 
higher virial densities (see our discussions in § [5] and § [3| ■ 
In §§ 1 4161 we quantify these effects and test the expected 
scalings. We verify that the analytic technique of Bullock et 
al. (2001, B01) for predicting halo concentrations works well 



when applied directly to constant w cosmologies, implying 
that it is fairly generally applicable. 

A related issue is the effect of dark energy on the halo 
mass function. Although it is not directly observable, the- 
oretical insight into the expected halo mass function is at- 
tainable through current N-body simulations. The effects 
of dark energy on halo mass functions have been investi- 
gated by Bartelmann, Perota, & Baccigalupi (2003), Linder 
& Jenkins (2003), Klypin et al. (2003), Maccio et al. (2003), 
and Lokas, Bode, & Hoffman (2003) and the accurate pre- 
diction of halo mass functions as well as accretion histories 
and density structures over a large range of redshifts is nec- 
essary in order to take full advantage of the ability of up- 
coming cluster surveys, such as the Sunyaev-Zeldovich Array 
(http://astro.uchicago.edu/sza/ ), to constrain the dark en- 
ergy equation of state (see Carlstrom, Holder, & Reese 2001 
for a review). There appears to be general agreement that 
halo mass functions can be approximated accurately by the 
"universal" mass function of Jenkins et al. (2001, hereafter 
J01) even in models with dark energy. In § 15.11 we confirm 
the results of previous studies of quintessence cosmologies 
with w > —1, and extend this agreement to models with 
w < -1. 

Klypin et al. (2003) and Dolag et al. (2003) have 
performed previous numerical studies of CDM haloes in 
quintessence cosmologies. Where our study overlaps with 
these, our results are generally in agreement. This study ex- 
tends, complements, and improves upon previous studies in 
several ways. We have investigated the effects of normal- 
izing the power spectrum using clusters vs. CMB, which 
yield quite different results. We have assembled large cata- 
logues of haloes, with masses and concentrations for more 
than 1600 haloes with more than 100 particles each in each 
simulation. With this large number of haloes we are able to 
measure the distribution of concentrations at fixed mass and 
test the predictions of the B01 model for both the mean re- 
lation between concentration and mass and the halo-to-halo 
scatter. We also extend the range of quintessence parameter 
space probed by simulations by exploring two models with 
w < -1. 

The format of this paper is as follows. After a brief 
review of the basics of structure formation in models with 
dark energy in §|21 we go on in §|3]to describe modifications 
to the B01 model for halo concentrations in quintessence 
cosmologies. In § 2] we describe the set of numerical sim- 
ulations that we performed. In § [5J we present the mass 
functions and concentrations of CDM haloes in these simu- 
lations and compare them to the results of the analytic B01 
model. In § |SJ we briefly discuss the central density prob- 
lem of CDM in light of our results. In § |7| we present a 
summary of our findings, highlight improvements and dif- 
ferences to previous studies, and discuss the implications of 
this work. Throughout this paper, we assume a flat universe 
with present matter density relative to critical of Qm = 0.3, 
and Hubble parameter h — 0.7. The quintessence energy 
density is then Qq = 1 — Qm = 0.7. We refer to a model 
with a cosmological constant (w = — 1) as ACDM. We use 
the terms "quintessence" and "QCDM" to refer to all mod- 
els which have w ^ — 1. 
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2 STRUCTURE FORMATION IN 
QUINTESSENCE COSMOLOGIES 

In this Section, we present a brief overview of the growth of 
structure in quintessence cosmologies. In H2.ll we introduce 
the basic results of cosmological perturbation theory with 
quintessence and discuss our computations of the matter 
power spectrum. We briefly discuss the normalization of the 
power spectrum in >I2,2I A convenient way to define the mass 
and radius of a dark matter halo is through a mean over- 
density A v i r , relative to the background density (see H3.1H . 
The idea is to choose A v j r so as to delineate the boundary 
between virialized and in-falling material. The equivalent 
linear overdensity at collapse 8 C , is a quantity that is used 
to delineate the mass scale of objects that are forming at 
a particular redshift. Often, these overdensity criteria are 
chosen with reference to the evolution of a spherical tophat 
overdensity. In H2.3I we discuss the spherical tophat model 
in quintessence cosmologies. 

We explore both conventional quintessence models with 
w > — 1 and, adopting an empirical approach, we pursue 
models with w < — 1. Requiring the kinetic energy of the 
quintessence to be positive imposes the condition w > — 1 
on the quintessence equation of state. When w > — 1 and 
constant, all computations can be performed knowing only 
the value of w by assuming that the scalar field Lagrangian 
has a canonical kinetic term. To explore the parameter space 
w < — 1 self-consistently, we must choose a particular model 
for scalar field dynamics and one simple possibility is a 
Lagrangian with a non-canonical kinetic term that differs 
from the canonical case by a negative sign. Explicitly, we 
adopt the "phantom energy" Lagrangian for the field, cj>: 
C = -d M 0<9 M 0- V{4>) (Caldwell 2002; Carroll, Hoffmann, & 
Trodden 2003; Cline, Jeon, & Moore 2003). Having made 
this choice, we can express the derivatives of the poten- 
tial completely in terms of w (Dave, Caldwell, & Steinhardt 
2002). 

2.1 Cosmological Perturbations 

Here we consider linear perturbations to the CDM density 
field, S(x) = Sp(x)/pb, with ph the mean density of dark 
matter in the universe. The linear evolution of 8(x) follows 
from solving the linearized Einstein-Boltzmann equations 
(e.g., Ma & Bertschinger 1995). We perform our calcula- 
tions in the synchronous gauge with a modified version of 
the publicly- available Einstein-Boltzmann code CMBfast by 
Seljak & Zaldarriaga (1996). We assume a scale-invariant 
spectrum of adiabatic primordial density fluctuations and a 
baryon density Q,bIi 2 = 0.02. 

The linear equation for fluctuations in the quintessence 
field in Fourier space is 

S<t> + 3HS<t> + (k 2 /a 2 ± V,^ )5<t> = S k [(l + w Q )p 4> ] 1/2 . / (1) 

where V,^ = d 2 V/d<j) 2 , is the mean energy density in the 
quintessence field, and the plus (minus) signs correspond to 
fields with positive (negative) kinetic energy. Here the dots 
represent derivatives with respect to cosmological proper 
time t, and Sk is the linear CDM fluctuation in Fourier space. 
The power spectrum of linear density perturbations is de- 
fined as P(k) = (\5k\ 2 )- The CDM transfer functions for our 
scale-invariant primordial spectrum are defined from 
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Figure 1. The ratio of z = quintessence transfer functions to 
ACDM for different models. The transfer functions are similar for 
wavenumbers k > 0.01 h Mpc" 1 . The different equation of state 
parameters are shown in the legend. 

P(k,z)=A Q kT 2 (k,z)^ y (2) 

where Aq is the normalization and D(z) = S(z)/S(0) is the 
linear growth factor. On scales larger than the Compton 
wavenumber, &q ~ y/V,^, 5(f> can grow and source the evo- 
lution of 5k- Contrarily, on small scales k > &q, perturba- 
tions in the quintessence field decay, so the QCDM transfer 
functions have the same form as those of standard ACDM, 
while on scales k < fcp, the transfer functions reflect the 
clustering of </> (Ma, Caldwell, Bode, & Wang 1999). Fig- 
ure shows the ratio of the CDM transfer functions with 
quintessence, Tq , to the transfer function in the correspond- 
ing ACDM model (w = -1), Ta, at z = 0. For k < k Q and 
w > — 1, perturbations in the quintessence field, 5p$ — 38p^,, 
source the growth of S^- For w < — 1, the source term from 
the quintessence field changes sign, resulting in a relative 
decrease in St- 

2.2 Power Spectrum Normalization 

One of the most important parameters determining the scale 
radii of halo density profiles is the normalization of the mat- 
ter power spectrum on small scales. In our N-body exper- 
iments, we choose to normalize each w model with similar 
values of erg, set by cluster abundance estimates. We do so 
because the abundance of clusters is a more direct probe of 
the amount of power on the scales that are relevant to galaxy 
formation than the CMB anisotropy measurements and in 
order to isolate the differences that arise because of varia- 
tions in the expansion rate in models with different values 
of w. 

We normalize the power spectrum to the abundance 
of massive clusters using x-ray flux and temperature mea- 
surements from the cores of clusters of galaxies, and the 
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corresponding conversion to cluster mass from the mass- 
temperature (M-T) relation. With all other cosmological 
parameters fixed, a comparison of the observed mass func- 
tion to the predicted mass functions from N-body simula- 
tions (J01) determines the normalization parameter erg. The 
quintessence field is smooth on scales much smaller than the 
horizon size, thus the values of erg derived from z = cluster 
measurements are rather insensitive to w. 

The largest source of error in the determination of erg 
from clusters is the uncertainty in the normalization of the 
M-T relation. Simulations consistently predict an M-T re- 
lation that is a factor of ~ 2 greater than the one observed 
from x-ray temperature data (e.g., Seljak 2002 and refer- 
ences therein). For ACDM and Qm = 0.3, our current un- 
derstanding of the M-T relation places erg very broadly in 
the range 0.65 <cr 8 < 1.1 (Pierpaoli, Scott, & White 2001; 
Reiprich & Bohringer 2002; Pierpaoli et al. 2003). Using the 
HIFLUGCS cluster mass function (Reiprich & Bohringer 
2002), we find for fijn = 0.3, erg ~ 0.74, nearly indepen- 
dently of w (Kuhlen, Keeton, & Madau 2004). We normalize 
the power spectra to as = 0.742, 0.740, 0.738, 0.736, 0.734 for 
w = —1.5, —1.25, —1.0, —0.75, —0.50 respectively. In princi- 
ple, the value of erg as determined from clusters is degenerate 
with f^M (e.g. Schuecker et al. 2003), indeed the global best 
fit to the HILFUGUS sample is Q M ^ 0.12, a 8 ~ 0.96 in a 
ACDM cosmology; however, we adopt values of erg that best 
fit the data given our choice of Qm = 0.3. 

While we have chosen to normalize our numerical sim- 
ulations using the cluster abundance, this normalization is 
fairly uncertain. In the interest of completeness, we remark 
that not all of these model normalizations are consistent 
with n = 1 normalization to CMB anisotropy (even modulo 
uncertainties in the reionization epoch). In the simplest case, 
with all other cosmological parameters held fixed, the value 
of w affects the CMB-derived erg normalization primarily 
through the late-time integrated Sachs- Wolfe (ISW) effect 
(see Hu & Sugyiyama 1995). For universes with w > — 1, 
the quintessence energy density becomes comparable to that 
of CDM at earlier epochs relative to ACDM, resulting in 
greater variation in the gravitational potentials along lines 
of sight to the surface of last scattering. For lower values of 
w, the ISW effect is not as prominent because quintessence 
becomes dynamically important only at more and more re- 
cent epochs. In Figure [5] we show the value of erg implied 
by the CMB normalization as a function of the equation of 
state parameter w. As we stated above, the decrease in erg 
is due to the increased importance of the ISW effect as w 
increases. 

Care must be taken when setting the CMB normaliza- 
tion. First, we note that when we performed our N-body 
experiments we fixed the parameter Dm — 0.3 for all models 
rather than allowing Qm to vary along the SIm-m degeneracy 
in the angular diameter distance. Our most extreme models 
are currently disfavored by the present observational data 
(e.g., Tegmark et al. 2003; Knop et al. 2003), but as our in- 
tent is to address the effect of quintessence on structure and 
halo formation, we do not consider this to be a serious defi- 
ciency. In principle, the WMAP (Spergel et al. 2003) result 
of high optical depth to the last scattering surface r, has 
made the determination of erg from the CMB less robust, 
as the scattering off of free electrons damps anisotropies 
on scales that are sub-horizon at the epoch of reionization, 
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Figure 2. The power spectrum normalization erg, implied by 
CMB anisotropy as a function of w. The solid line shows the 
values of trg that we obtain by assuming the optical depth to the 
last scattering surface r = 0. The dashed line shows the values 
of erg implied by adopting r = 0.17. The dotted line shows the 
values of erg that we infer from the abundance of massive x-ray 
clusters. 



thereby introducing a degeneracy between r and erg. The 
dashed line in Figure shows the CMB-normalized erg as a 
function of w implied by adopting an optical depth to the 
last scattering surface of r = 0.17, in line with the WMAP 
expectations. 1 We return to a discussion of the relative im- 
portance of erg and w in §5 16171 



2.3 The Spherical Collapse Approximation 

A common convention is to define the virial mass and radius 
of a dark matter halo by demanding that the mean density 
within the virial radius of the halo be a factor A v i r times 
larger than the background density, pb ■ Thus the virial mass 
and radius of a halo are related by 



M vil 



47V A p3 
— /AvirPb^tvir 



(3) 



In addition, the equivalent linear overdensity at collapse 
5 c (z) is often used to determine the mass scale that is typ- 
ically collapsing at a given epoch. Both of these quantities 
are usually estimated using the approximation of spherical 
tophat collapse (e.g., Lacey & Cole 1993). In this section, 
we summarize the previous results for spherical collapse in 
quintessence cosmologies (Mainini et al. 2003; Mota & van 



1 Note that a high optical depth to reionization appears to be dif- 
ficult to reconcile with low values of erg < 0.75 (e.g., Somerville, 
Bullock, & Livio 2003), though the CMB-derived r is very un- 
certain and is somewhat degenerate with the tilt of the power 
spectrum, so tilted models with low erg may still be viable. This 
issue will likely be settled with future CMB data and analysis. 
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Figure 3. Ratio of the growth factor to the SCDM growth factor 
(D scdm(i) = a) as a function of scalefactor for different values 
of w. The growth factor is normalized to unity today (D(a = 1) = 
1). Linetypes are shown in the legend. 



de Bruck 2004), and discuss the implications for w < — 1 
cosmologies. 

The evolution of linear overdensities on scales much 
smaller than those on which the quintessence field spatially 
clusters is 

3 



5 + 2H(a)S~ -H Q 



Ma 



(4) 



Solving this equation gives the growth factor for linear per- 
turbations. Figure |21 shows the linear growth factor in five 
different quintessence models, normalized to the growth fac- 
tor in an $1m = 1, SCDM cosmology in which D(a) oc a. 
Models with w > — 1 have been well-studied, with such 
models showing more relative growth at higher redshifts. 
For models with w < —1, figure |3] shows how the trend to- 
wards more relative growth at lower redshifts continues for 
w < -1. 

To detemine the non-linear growth of an object which 
has decoupled from the expanding universe and virialized, 
we follow Wang & Steinhardt (1998) and Weinberg & 
Kamionkowski (2003). We use the spherical collapse ap- 
proximation to determine the non-linear overdensity of a 
halo, Avir(z), as a function of w. In this model, an over- 
density defined by A v i r (z) collapses at the time t co n, when 
the radius of the overdense region approaches zero. How- 
ever, the actual, final radius of the collapsed object is finite 
and can be computed using the virial theorem. We compute 
the equivalent linear overdensity at collapse by evolving the 
linearized equation of motion, Eq. g|, until the time t co ii, de- 
termined from the non-linear evolution of the overdensity. In 
the ACDM cosmology, the well-known results for the linear 
and non-linear overdensities at collapse are S c (z = 0) ~ 1.67 
and A v i r ~ 337, and vary with redshift. 

Figure[I]shows our results for the equivalent linear over- 
density at collapse 8 c (z), and the non-linear overdensity at 
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Figure 4. The non-linear and linear overdensities at collapse in 
quintessence cosmologies. In panel (a), we exhibit the non- linear 
overdensity at collapse A v ; r (z) for 5 different quintessence models 
that we explore in this paper ((1m = 0.3). The different values 
of the quintessence equation of state parameter w, are shown in 
the legend. The increase in A v - lr (z) for larger w results primarily 
from later collapse redshifts. In panel (b), we show the equivalent 
linear overdensity at collapse S c (z), as a function of redshift for 
the same quintessence models. 



collapse Avir(z), as a function of redshift in quintessence 
cosmologies. The trend in 5 c (z) reflects the fact that over- 
densities grow more slowly in higher-ui models. We can un- 
derstand the behavior of A v i T (w, z) by noting that overden- 
sities in models with larger w take longer to collapse at late 
times. Consider, for example, two haloes, one in a w = — 1 
cosmology, and another in a m > —1 cosmology, both of 
which are just virializing at some redshift z. Very roughly 
speaking, the w > — 1 halo will have had a turn-around time 
at a higher redshift than the w = — 1 halo, and its virial den- 
sity will be higher to reflect this. Conversely w < — 1 models 
will have more recent turn-around times and smaller virial 
overdensities. We find that an accurate fitting function, in- 
cluding the regime w < — 1, can be obtained from a slight 
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modification to the formula already proposed by Weinberg 
& Kamionkowski (2003) for w > — 1, 



A vir (z) = 18tt 2 [l + ae b (z) 



(5) 



where @(z) 

|0.234 



0.432 



c M y^./ — 1, and with a 
2.001(|«| u -"* - 1) and b 
We find this formula to be accurate to better than 2% for 
0.1 < fi M < 1 and -0.5 < w < -1.5. 



0.929 - 0.222(H ' 727 - 1) 



3 ANALYTIC MODEL FOR HALO 
CONCENTRATIONS 

3.1 Main Ingredients 

It is commonly agreed that the spherically-averaged density 
profiles of dark matter haloes can be described fairly well 
by a generalized NFW profile on scales that are resolved in 
state-of-the-art N-body simulations: 



Phalo(f) 



Ps 



(r/r s ) a (l - r/r s ) 3 



(6) 



where a describes the slope of the inner density profile at 
r < r s . The value of a that most closely represents the re- 
sults of N-body simulations is still debated, with acceptable 
values between —0.7 and —1.5. An additional complexity 
is that recent studies indicate that haloes exhibit a range 
of inner slopes (Klypin et al. 2001; Tasitsiomi et al. 2003; 
Navarro et al. 2003). In the following, we adopt a = — 1, cor- 
responding to the standard NFW profile. In this paper, we 
are concerned with the concentration parameter c v i r , which 
is quite insensitive to the exact value of a (see below). 

The two parameters of the NFW profile are r a and p a , 
with r s the radius at which the logarithmic density slope 
becomes equal to —2. The concentration of the halo is de- 
fined as the ratio of its virial radius to the scale radius of 
the NFW profile, 



r B 



(7) 



With these definitions, the halo virial mass is related to the 
NFW parameters by 



A/vir = 4np s r 3 3 



ln(l + c vir ) - 



1 + Cvir 



(8) 



so the halo density profile is completely determined by M v ir 
and Cvir. 

Numerical simulations have revealed a correlation be- 
tween Mvir and c v i r , with halo concentrations log- normally 
distributed around the median relation. Several simple mod- 
els have been developed to explain this correlation (NFW; 
B01; Eke, Navarro, & Steinmetz 2001). Here we focus on the 
B01 model and test the accuracy with which it predicts the 
observed relation between halo mass and concentration in 
simulations with w ^ — 1. 

The B01 model assumes that a halo's central density 
and concentration are set by the density of the universe at a 
characteristic formation epoch. This formation epoch quali- 
tatively tracks the characteristic collapse epoch for the halo 
subunits. It is defined as the time when the linear rms den- 
sity fluctuations at a scale corresponding to a fraction F of 
M v ir is equal to the linear collapse overdensity, 8 C : 



Given the halo collapse epoch, the halo concentration 
is set via 

a 



Cvir (Mvir, a) = K x 



a c (M v 



(10) 



A closely related study by Wechsler et al. (2002, here- 
after W02) showed that the B01 "collapse epoch" seems to 
correspond closely to the epoch when the mass accretion rate 
of the halo d ln(M v ir)/d£, is large compared to the rate of 
cosmic expansion. W02 found that if one defines the end of 
this rapid collapse phase to be when dln(M v i r )/dln(a) < 2, 
it corresponds closely to the "formation epoch" in Eq. @ 
above. After this epoch of rapid mass accretion ends, the 
halo mass and virial radius continue to grow via compara- 
bly minor mergers and diffuse mass accretion. These rela- 
tively minor mergers do not affect the inner regions of the 
halo (r < r„) significantly and so R v i T grows, but r s remains 
approximately constant, leading to an increase in concentra- 
tion as the halo evolves (W02). 

The B01 model has two parameters, F and K, that have 
to be determined by calibrating them to numerical simula- 
tions. B01 analysed two ACDM simulations with erg = 1.0 
and different resolutions and box sizes. Using a halo finder 
based on the Bound Density Maxima (BDM) algorithm 
(Klypin & Holtzman 1997), they assembled a catalogue of 
several thousand haloes for each simulation, and fit NFW 
profiles to each of them. B01 found that their model was 
able to satisfactorily reproduce the mean relation between 
halo mass and concentration, and the redshift dependence 
of the Cvir-Mvir relation with F = 0.01 and K = 4.0. 
B01 and W02 determined that the scatter in concentration 
at fixed mass is well-described by a log-normal distribu- 
tion with oiogc = 0.14 dex. Other numerical studies have 
found a somewhat smaller scatter, with Jing (2000) report- 
ing aiogc = 0.9 — 0.11 dex, and Jing & Suto (2002) finding 



0.13 dex. 



a(FM viT ,a c ) = S c (a c ). 



(9) 



The B01 model with F = 0.01 and K = 4.0 has since 
proven successful in reproducing concentrations of ACDM 
haloes over more than six orders of magnitude in halo mass 
from A/vir ^ 10 7 M Q to M vir ^ 10 13 M (e.g., Colin et al. 
2003; Hayashi et al. 2003). However, as discussed in B01, the 
model ceases to make physical sense for halo masses large 
enough that FM V1I begins to approach the typical collapse 
mass at z — 0. This is because linear fluctuations stop grow- 
ing at late times in ACDM, and with the simplified defini- 
tion of collapse time discussed above, very large haloes never 
collapse. Consequently, the B01 model with F = 0.01 under- 
predicts halo concentrations for systems more massive than 
a few x 10 14 M (e.g., Hayashi et al. 2003; Dolag et al. 2003). 
As a matter of pragmatism, this can be remedied with a sim- 
ple change of parameters. With F = 0.001 and K = 3.0 the 
model works adequately for all masses, though it becomes 
somewhat less attractive because of the small value of F, for 
which the collapse epoch no longer corresponds directly to 
a c defined by W02. 



3.2 The Analytic Model in Quintessence 
Cosmologies 

As shown in § |21 linear overdensities have greater relative 
growth at higher redshift as w increases. We then expect, 
given an overdensity on a mass scale M v i r , that this mass 
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scale will collapse at higher redshift as we increase w. As 
the halo concentration reflects the density of the universe at 
the time of rapid collapse, we expect this change in average 
formation time to translate directly into a change in average 
concentration. The differences in to 7^ —1 models should be 
confined to changes in the rapid-collapse epoch a c for haloes 
of a given mass. 

From Equation JUJ, a c for a halo of mass M v i T is de- 
termined by S c , D(a), and erg. The changes in S c with w 
(Figure 01 have a very small effect on a c . More relevant are 
the changes in a(M,a) and the linear growth rate (Fig. |3J. 
For a fixed value of erg all the ui-dependence of a c (M v j r ) will 
be captured by D(a), with a(M,a) reaching the S c collapse 
threshold at earlier times as w increases. We therefore expect 
concentrations to increase as w increases. 

Note that an alternative definition of a halo's radius, 
-R200, the radius at which the mean halo density is equal to 
200/Ocrit has frequently been used in the past. This results 
in an alternative definition of concentration: C200 = -R200 /r s 
(e.g., NFW96). Of course, given a set of cosmological pa- 
rameters, the model described above (and similar models) 
can be used to predict equivalent relations between C200 and 
M200 • For a fixed cosmology, the predicted relation between 
C200 and M200 will look quite similar to the c v i T -M v i T relation, 
with an offset that varies slowly as a function of concentra- 
tion and accounts for the differences in values of the outer 
halo radius. 

Because the simulations of B01 focused only on one cos- 
mology, it was impossible to tell whether agreement with the 
simulations and the proposed B01 model was sensitive to the 
choice of defining halo concentration relative to i? v j r instead 
of i?200- That is, one could have equally well proposed a 
different model based on C200: 

C2oo(M 20 o, a) — K200 — 777 — c- (H) 
a c (M 2 oo) 



The simulation results of B01 could have been reproduced 
using this model, simply by setting K200 equal to a slightly 
smaller value than the original K = 4.0. 

Consider now the current case, where we compare sim- 
ulation results based on cosmologies with different w's and 
hence different A v ir values. In these simulations, the ratios 
of -Rvir/^200 (and Cvi r /c2oo) f° r fixed-mass haloes will de- 
pend on the value A v i T (w). Therefore, it is impossible for a 
model based on C200 with fixed K200 to do equally well as 
the original B01 model based on c v j r with fixed K for all val- 
ues of w. More physically, the original B01 model described 
above implicitly assumes that the haloes have virialized at 
the appropriate virial density and predicts that, in addi- 
tion, the halo collapse redshift acts to set the ratio of the 
virial density (oc R~?) to the central density (cx r~ 3 ). A 
model based on C200 with fixed K200 would instead assume 
that all of the changes in halo density arise solely because 
of changes in a c . As we demonstrate below, the virial as- 
sumption seems to capture better our simulation results. It 
seems therefore, that there are two physical processes that 
set halo densities: one process is related to the global pro- 
cess of halo virialization and the other may be related to an 
earlier, rapid-collapse epoch. 



4 NUMERICAL SIMULATIONS 

In this section, we describe our numerical simulations. In 
§ 14.11 we detail the numerical and cosmological parameters 
that were used. In § 14.31 we describe the methods we use 
to locate haloes and to fit NFW profiles to their density 
profiles. 



4.1 Simulations and Parameters 

We use GADGET version 1.1, a publicly-available and well- 
tested iV-body code (Springel et al. 2001). Gravity between 
particles is solved using a hierarchical tree algorithm in co- 
moving coordinates, and both the force calculations and the 
time-stepping are performed in a fully adaptive way. Us- 
ing the parallel version, we have run the code on either 96 
375MHz IBM Power3 processors of NERSC's Seaborg or on 
64 1.4GHz Athlon processors of UpsAnd, a 264-processor Be- 
owulf cluster at The University of California at Santa Cruz. 
We made necessary alterations to the expansion rate of the 
universe for GADGET to account for quintessence cosmolo- 
gies with w =fc — 1. 

Power et al. (2003) have performed a detailed conver- 
gence study of a high resolution cluster simulation using 
GADGET, and although we simulate a much larger cos- 
mological volume we have followed their recommendations 
for a number of GADGET'S parameters. In particul ar we 
have chosen an adaptive timestep equal to Ati = r\ a t \ftifai, 
where e; and Oi are the gravitational softening and accelera- 
tion experienced by the i th particle in the simulation, and rj at 
is a dimensionless constant. Power et al. (2003) recommend 
setting r\ae = 0.2; this choice of adaptive timestep minimizes 
undesirable effects due to particle discreteness and hard 
scatterings, while at the same time allowing for convergence 
at minimal computational expense. In GADGET, gravita- 
tional softening is performed using a cubic spline (Springel 
et al. 2001), for which the potential becomes exactly Newto- 
nian at r — 2.8 e, where e is the softening length. Generally 
our simulations were run with a co-moving softening length 
of e = 2.5/i _1 kpc, although we have run a few cases with e 
as low as Ih" 1 kpc. 

Our cosmological background model is fixed by Qm — 
0.3, S1q = 0.7, h = 0.7, and n = 1.0 for all values of w. 
In normalizing erg on the scale of galaxy clusters, the initial 
power spectra are nearly unaffected by quintessence. How- 
ever, when normalizing to the scales probed by the CMB, 
the initial power spectra are altered by the inclusion of 
quintessence (section 12. 21 . As discussed below, we normal- 
ize our simulations such that erg ~ 0.74, thus the effect of 
w ^ — 1 is due almost exclusively to the expansion rate. 

All of our simulations were run with 256 3 particles in 
boxes with sides of length 60/i~ 1 Mpc. Qm = 0.3 implies a 
mass per particle of M v = 1.1 x 10 9 /i _1 Mq. For the analysis 
of halo concentrations we used only haloes with more than 
100 particles (see *I5.2I| . This corresponds to a minimum halo 
mass of M^ n ~ 1.1 x lO 11 /! -1 Mq, for which the B01 model 
predicts a median concentration of 13.5 (for erg = 0.74). This 
translates into an NFW scale radius of r™ m ~ 10/i -1 kpc. 
More massive haloes will have larger scale radii, and because 
even this minimum scale radius is almost three times larger 
than our softening length, we should be able to determine 
accurate concentrations from the haloes in our simulations. 
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Table 1. Simulation Parameters and Power Spectrum Normal- 
izations. 



Model w 


e/h 1 kpc 


°"8,nom 


0"8,off a 


-0.50 


2.5 


0.742 


0.799 


-0.75 


2.5 


0.740 


0.775 


-1.00 


2.5 


0.738 


0.716 


-1.00 


1.0 


1.000 


0.972 


-1.25 


2.5 


0.736 


0.716 


-1.50 


2.5 


0.734 


0.714 



All other parameters are fixed at the same value for all simu- 
lations. The number of particles is N p = 256 3 , the box size is 
Lbox = 60/i — ^Mpc, and the initial redshift is Zi = 50. For all sim- 
ulations, the remaining cosmological parameters are Qm = 0.3, 
Uq = 0.7, h = 0.7, and n = 1.0. 

a the difference between o^nom and Grg, e ff is explained in -14.21 

Table summarizes the parameters used in our simula- 
tions. 



100 




Figure 5. A histogram of the effective values of as determined 
from 1000 Monte-Carlo realizations of the initial conditions with 
an input powcrspcctrum normalized to erg = 0.74. The distribu- 
tion is Gaussian with a standard deviation of 0.049, or roughly 
7%. 



4.2 Initial Conditions 

Setting initial conditions for our simulations requires fixing 
the z = power spectrum normalization, which we param- 
eterize by erg. As discussed above f ij2.2l . our current under- 
standing of present day clusters of galaxies makes erg still un- 
certain by ~ 20-30% (erg ~ 0.70-1.10 for fi M = 0.3). Halo 
concentrations in cosmological N-body simulations depend 
sensitively on erg, because the amount of small-scale power 
directly affects typical halo formation times, especially for 
high mass haloes. For example, at M v i r = 10 14 h" 1 Mq the 
B01 model predicts median concentrations of c v j r = 6.8 
and Cvir = 5.3 for erg = 0.90 and erg = 0.74, respectively. 
We initialize our simulations with the values of erg deter- 
mined by Kuhlen et al. (2004) from the abundance of clus- 
ters in the HIFLUGCS sample of local clusters (Reiprich 
& Bohringer 2002): erg = 0.742,0.740,0.738,0.736,0.734 
for w = -1.50, -1.25, -1.00, -0.75, -0.50, respectively. We 
construct initial conditions for each w with the routines of 
the publicly-available PM code (Klypin & Holtzman 1997). 

The process of initializing particle positions and veloci- 
ties based on the linear power spectrum is subject to statisti- 
cal fluctuations. The largest modes of the system (A ~ I/box) 
are sampled only twice per dimension and are thus sen- 
sitive to deviations caused by small number statistics. As 
8/i _1 Mpc is close to Lbox, this can lead to noticeable differ- 
ences between the "nominal" value of crg >nom , used to con- 
struct the initial conditions and an "effective" value of crg^ff, 
determined from the actual particle positions at the initial 
redshift. In order to quantify this difference, we determine 
°"8,off by direct numerical integration of the initial A-body 
power spectrum: 



2 1 D l (z 

°"8,cff — 



0) 



2tt 2 D 2 (z = 50) 



,{k)W{kRg,y&k, (12) 



where P n um(&) is the power spectrum derived from the N- 
body initial conditions and W(x) is the spherical tophat 
window function given by W(x) = 3/x 2 {sm.x/x — cos a;), 
evaluated at kRg, with i?g = 8/iT 1 Mpc. 

Pnum(k) only extends to fc m i n w 0.1h~ Mpc -1 , which 



is not low enough to allow the integral in Equation 1121 to 
converge. We estimated the portion of the integral below 
fcmin by integrating the smooth analytical power spectrum, 
and applied this correction to get crg ie ff ■ We found that for 
w < —1 CTg jC ff ~ 0.715, which is ~ 3% lower than crg, n om. 
However, the variation of erg in a box of this size due to 
cosmic variance should be roughly ~ 5%, so this difference 
is not surprising. To demonstrate this explicitly, we have 
constructed 1000 realizations of the initial conditions and 
computed values of erg. We infer from these realizations that 
measured values of the effective erg are distributed with a 
standard deviation of ~ 5%. The resulting distribution is 
shown in Fig.[£] Therefore, the difference between erg, e ff and 
erg,nom for the three w < — 1 cases is not surprising and is 
consistent with cosmic variance. We find similar deviations 
of erg^eff from erg inom for all of the values of w that we simu- 
late. 



4.3 Halo Finders 

We use two different halo finding algorithms to locate 
the haloes in our simulations, depending upon the quan- 
tities we probe with our simulation data. In § 15.11 
we compare the mass functions of our simulated haloes 
with the J01 "universal" mass function. J01 used the 
Friends-Of-Friends (FoF) algorithm (Davis et al. 1985) 
to identify simulated haloes. In order to make a di- 
rect comparison to the J01 mass functions, we have 
employed a University of Washington FoF halo finder 
( http : //www-hpec . astro .Washington. edu/tools/f of . html ). As 
in J01 we set the linking length to 0.2 times the mean inter- 
particle separation for all models. 

To directly compare our halo concentrations to the B01 
model (§ I5.2|l . we use an updated version of the halo finder 
that B01 and W02 employed in order to identify haloes. This 
halo finder is based on the BDM algorithm (Klypin & Holtz- 
man 1997) and iteratively removes particles that are not 
bound to the halo in question. Upon identifying haloes, we 
fit NFW profiles to each halo and determine c v i r - For more 
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detail, we refer the reader to Appendix B of B01 and Ap- 
pendix A of W02. We include in our catalogues only haloes 
with more than 100 particles, the same cut-off used by W02. 

We have checked that both halo finding algorithms 
agree with each other, within our expectations, by compar- 
ing mass functions. The systematic differences in total mass 
between haloes defined in terms of the cosmology-dependent 
virial overdensity A v i r (w) (as in the BDM finder) and those 
based on the fixed FoF linking length translate into only 
minor differences in mass functions. At high redshift, the 
two mass functions actually converge because A vlr (z) ap- 
proaches 178 and a linking length of 0.2 times the mean 
inter-particle separation roughly corresponds to a mean halo 
density of ~ 180 — 200 times the background density. Below 
z ~ 2 — 2.5, the FoF and BDM mass functions agree well. 
At higher redshifts the BDM-based finder becomes increas- 
ingly incomplete. A consequence of the low value of erg in 
our simulations is that halo formation occurs more recently. 
Frequent merger events during the rapid mass growth phase 
of halo formation disrupt any spherical symmetry in the halo 
density profile. These haloes will not be well described by 
the NFW formula. In our implementation of the BDM halo 
finder, haloes with very bad fits to the NFW profile are re- 
jected and not included in the catalogues. This is the major 
source of incompleteness at z > 2.5. Note that the net effect 
of this incompleteness is to underestimate the number of low 
concentration haloes at high redshift. 



5 RESULTS 

5.1 Mass Functions 

Several recent numerical studies have demonstrated that the 
J01 formula for halo mass functions may be considered "uni- 
versal" as it accurately describes halo counts as a function 
of mass in iV-body simulations of various cosmologies, in- 
cluding models containing dark energy with w ^ — 1 and a 
time-varying w (Linder & Jenkins 2003; Klypin et al. 2003; 
Maccio et al. 2003; Lokas, Bode, & Hoffman 2003). We con- 
firm and further extend this conclusion by presenting halo 
mass functions from our simulations at different redshifts. 
In Figure |S] we show the mass functions of our FoF haloes 
in each cosmology and at a variety of redshifts from z = to 
z = 3. It is apparent that, in all panels, the mass functions 
and the redshift evolution of the mass functions for each 
w model are in excellent agreement with the J01 forumula. 
Thus we confirm that the J01 formula is a good approxima- 
tion at all redshifts for w > — 1 and we extend the range of 
validity of the J01 relation to include quintessence models 
with w < — 1. 



5.2 Concentrations 

In H3.2I we described the manner in which quintessence mod- 
ifies the predictions of the analytic B01 model for halo con- 
centration as a function of mass. We have assembled cata- 
logues consisting of more than ~ 1600 haloes for each of our 
quintessence iV-body simulations. The density profiles of ev- 
ery halo have been fit to an NFW profile, yielding a best-fit 
c v ir for each object. Some of these fits produced concentra- 



tions smaller than one. We have excluded these haloes from 
our subsequent analysis. 

The resulting c v ir(M v i r ) relations are plotted in Fig- 
ure [7| We find that our simulations produce haloes with 
slightly lower concentrations than expected from previous 
simulation results (e.g., B01, Colin et al. 2003) and the an- 
alytic model proposed by B01 with F = 0.01 and K — 4.0. 
However, we find that this difference can be described well 
by a constant offset. For example, keeping F = 0.01 and 
lowering the proportionality constant K to K — 3.5 [see 
Eq. 1101 in the B01 model matches our data quite well for 
all of the w models we explore. We discuss this overall off- 
set further in § 15.31 The c v i r -M v i r relation flattens out be- 
low M v ir ~ 6 x 10 u ft -1 M . We attribute this to the lower 
number of particles in these haloes, making them more sus- 
ceptible to relaxation effects which tend to cause the cen- 
tral regions of haloes to be more diffuse and lead to lower 
concentrations. It is thus unlikely that this flattening rep- 
resents any physical effect (compare to the results of Colin 
et al. 2003), and we have neglected the lowest mass bin in 
determining the best-fit value of K. 

As mentioned above ( H3.2H . considering several cosmo- 
logical models with different virial overdensities A v i r allows 
us to distinguish between analytic prescriptions based on 
definitions of halo concentration in terms of -R200 , in which 
the proportionality constant K200 is independent of cosmol- 
ogy (Eq. lilt , and those based on the virial radius i? v ir, 
in which K v [ r is cosmology- independent (Eq. 1101 . To test 
this, we re-analysed the five z = iV-body outputs us- 
ing the BDM halo finder, but setting A 

vir — pvii 

Pvir/pcrit £l M = 200 f2 M ~ 667, effectively yielding a 
relation between C200 and M2oo- Matching these relations 
to the model described by Eq. 1111 we determined best- 
fitting values of #200 = (3.76,3.44,3.32,3.16,3.16) for w = 
(-0.50, -0.75, -1.00, -1.25, -1.50), respecively. This range 
in -K200 is not consistent with one cosmology-independent 
value of -K200 ■ The results of this analysis suggest that mod- 
els similar to the B01 model, in which the halo concentration 
is defined in terms of _R v i r and A v i r , are more readily general- 
izable to alternative cosmologies cLS Cvir IS related to a/a c via 
a cosmology-independent constant of proportionality, K vlr . 
Put another way, defining the radius of a halo, and thus its 
concentration, using a fixed overdensity criterion necessi- 
tates using a cosmology-dependent proportionality constant 
in Eq. 110H while the cosmology-dependent virial overden- 
sity definition seems to account for these differences, so that 
-KVir is independent of cosmology. 

As in previous studies (B01; Jing 2000; Jing & Suto 
2002), we also find that haloes of a given mass have a broad 
distribution of concentrations. To determine the inherent 
scatter in the c v ir-M V ir relation it is important to account 
for the artificial scatter introduced by uncertainties in the 
fit to an NFW profile and by the Poisson noise in each bin. 
Following the B01 analysis, we corrected for the former by 
determining 500 one-sided Gaussian deviates for each halo 
with a standard deviation equal to the error in the c v i r fit 
returned by the halo finder. The deviates are positive (nega- 
tive) if Cvir is less (greater) than the median in that bin. We 
then determined the 16 th and 84 th percentiles in log(c v ir) 
and subtract off the Poisson noise from each in quadrature. 
The resulting estimates of the intrinsic scatter are shown as 
the dashed lines in Figure |7| The scatter is consistent with 
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Figure 6. A comparison between the J01 analytic mass function (solid lines) and the mass functions derived from our five N-body 
simulations (shapes with error bars, see the legend in the lower right portion of the Figure). The haloes in our simulations were located 
using a FoF algorithm with a linking length equal to 0.2 times the mean particle separation as in J01. In each panel we plot mass 
functions at various redshifts, from top to bottom: z = 0.0,0.5, 1.0,2.0, and 3.0. The error bars represent the Poisson noise due to the 
finite number of haloes in each mass bin. 



being independent of w and M v i r , and we find that tak- 
ing the B01 proportionality constant to be K\ ovl — 2.28 and 
^high = 4.52 fits the lower and upper lines well. These values 
correspond to <Ti og c j ow = 0.18 dex and o"i og c,high = 0.11 dex. 
Although these are similar in magnitude to the scatter re- 
ported in previous studies, our distributions are skewed away 
from log-normal toward lower concentrations. We note that 
the skewness may likely be caused by the lower resolution of 
our simulations, which tends to result in lower concentration 
haloes. 



For a fixed mass the B01 model predicts that concentra- 
tion should decrease with redshift as 1/(1 + 2). The haloes 
in our simulation also satisfy this relation, as shown by Fig- 
ure[5| in which we plot the redshift dependence of concentra- 
tion for haloes of mass M v ir = 7 x lO 11 ^ 1 M©. This figure 
shows that the concentrations follow the c vlr cc (1 + z) -1 
relation that is embodied in the B01 analytic model. At 
redshifts greater than ~ 2.5, our catalogues of haloes in 
this mass bin with fitted NFW profiles become incomplete. 
This incompleteness preferentially affects low concentrations 
haloes, causing the c v i r (z) relation to flatten at high redshift. 
We do not believe this to be a physical effect, and trust our 
data points to z ~ 2.5. 



5.3 The Concentration Discrepancy 

We have attempted to understand the origin of the discrep- 
ancy between c v ir(M v j r ) derived from our simulations and 
those reported by B01 and summarized by the B01 model. 
We have re-analysed the same z = simulation data that 
was analysed previously by B01, and we were able to repro- 
duce their c v j r -M v i r relation and scatter. We conclude that 
the discrepancy that we observe is not due to any change in 
analysis procedures. 

Of course, the main difference between the study of B01 
and our work is the choice of simulation codes. Whereas 
we use the publicly- available, uniform-resolution code GAD- 
GET, B01 used the adaptive-refinement code ART. Un- 
doubtedly the effective resolution at the centres of haloes 
was higher in the B01 simulation than in ours. In order 
to shed further light on this matter, we have run one ad- 
ditional GADGET simulation designed to test the impor- 
tance of the effective force resolution. Compared to the five 
simulations discussed previously, this one has higher force 
resolution (e = 1.0ft -1 kpc) and as = 1.0. The resulting 
Cvir(Af v i r ) relation is shown in Figure^ 

Here, we again find that the GADGET concentrations 
are systematically lower than the ones found by B01 with 
ART by - 14%. Instead of K = 4.0 we find K = 3.44 
matches the GADGET halo concentrations. This is con- 
sistent with the value found for the five lower resolution 
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Figure 7. The concentration parameter c v i r , as a function of mass M v i r , for the five quintessence N-body simulations. The original 
(F = 0.01, K = 4.0) and our best-fitting (F = 0.01, K = 3.5) B01 models are over-plotted as thin dotted and solid lines, respectively. 
Error bars represent the Poisson noise due to the finite number of haloes per bin. The dashed lines are our estimates of the intrinsic 
scatter in the relation, obtained by removing the scatter due to errors in the fits of the NFW profiles, as well as Poisson noise. The 
lowest mass bin was excluded in our analysis and is shown here for completeness only. The total number of haloes in the remaining bins 
is shown in the upper right corner of each panel. The lower right panel shows the best-fitting B01 model for all five values of w. 



quintessence simulations described above. The difference in 
K between our simulations and the B01 simulation may be 
due to an inherent difference between the GADGET and 
ART iV-body codes. Whether this is due to the higher max- 
imum force resolution afforded by the adaptive refinement of 
ART, or another difference between the two codes remains 
unclear. 



Several recent analyses based on w = —1 simulations 
with higher resolution than our own (Hayashi et al. 2003; 
Colm et al. 2003; Tasitsiomi et al. 2003) also favor the B01 
model with K = 4.0. In light of these results we suggest that 
our GADGET simulations systematically under-predict halo 
concentrations by ~ 10 — 15%. However, when this offset is 
normalized out, the variation of c v i r (M v ir) with w scales as 
predicted, and we conclude that the B01 model is successful 
in this regard. Emboldened by this success, we use the model 
to explore the implications of various normalization choices 
and to compare expected halo densities with those inferred 
from galaxy rotation curves. In what follows, we assume K 
1.0 for the model normalization and we advocate this choice 
for the reasons outlined above. (However, setting K — 3.5 
would not qualitatively change the conclusions that follow.) 



6 HALO CONCENTRATIONS AND CENTRAL 
DENSITIES WITH w ^ -1 

In Figure 1101 we illustrate the degeneracy between w and 
erg in setting halo concentrations. Shown are the predictions 
of the B01 model (F = 0.01, K = 4.0) for c vir (z = 0) with 
a fixed normalization as = 0.74 for several values of w. As 
discussed in § 3.2, for a fixed normalization, concentrations 
increase as w increases because haloes collapse earlier. The 
right panel of the figure shows the corresponding predictions 
for c v i r with erg determined by normalizing to the CMB with 
t — (see § 2.3). Note that for the CMB normalization, the 
trend with w is in the opposite direction, with increasing 
w implying lower concentrations. As illustrated in Figure |5J 
higher w requires a lower normalization: as — 0.6, 0.9, and 
1.0 for w — —0.5, —1, and —1.5, respectively. The change in 
normalization based on the CMB dominates changes in the 
growth function that give rise to the behavior of c v i T (w) at 
fixed normalization. 

As discussed previously, the c v ir parameter is useful, but 
it is not a direct measure of physical density. We would like, 
therefore to convert our predicted c v j r relations into quan- 
tities that have a more direct physical interpretation, and 
are more amenable to comparison with observations. Alam, 
Bullock, & Weinberg (2002) proposed the central density pa- 
rameter as a means to quantify the physical density in the 
central regions of a galaxy: 
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Figure 8. The median concentrations of haloes in bins of mass M vir = 7 x lO 11 ^- 1 M Q ±0.05 dcx as a function of rcdshift. Error bars 
represent the Poisson noise due to the finite number of haloes per bin. c v \ T falls in proportion to 1/(1 + z) (solid lines), in agreement 
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Ay/2 is the mean overdensity within ry/2, the radius at 
which the galaxy rotation curve reaches half its maximum, 

Knax ■ 

The A v /2 parameter is advantageous for several rea- 
sons. First, it facilitates comparisons between theory and 
observation. Any predicted c v i r vs. M v i r relation can be eas- 
ily converted into a A v /2 vs. l/ max relation. Similarly, given 
an observed galaxy rotation curve, A v /2 can be determined 
without reference to any particular analytic density profile. 
It also has the useful characteristic that even if an observed 
rotation curve is rising at the last measured point, substitut- 
ing the highest (outer) most point on the rotation curve for 
Vmax in the formula above results in an upper limit on the 
true value of A v /2- Specifically, if V max is underestimated 
by a factor / v , and the density profile varies with radius as 
p(r) oc r~ a , this leads to an overestimate of A v /2 by a fac- 
tor of j- 2q /( 2_q ) Thus this is an overestimate so long as 
the density profile falls off with radius or is constant. Fur- 
thermore, if a > 2/3, the fractional overestimate of A V /2 is 
larger than the fractional underestimate of Vmax- 

We compare the B01 model predictions in terms of these 
quantities in Figure [TT1 for three representative quintessence 
cosmologies to the observational data of low-surface bright- 
ness and dwarf galaxies compiled by Zentner & Bullock 
(2002; the same set of observed galaxies are discussed in 
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halo concentration changes when we normalize based on the CMB: erg ~ 0.6 , 0.9, and 1.0 for w = —0.5,-1, and —1.5, respectively. 



Hayashi et al. 2003, who reach similar conclusions) from the 
observational work of Swaters (1999), de Blok, McGaugh, & 
Rubin (2001), and de Blok and Bosma (2002). The error bar 
in the right panel shows the theoretical lcr scatter in Ay/ 2 
expected at fixed Vmax due to the scatter in c v i r . 

While it is difficult to quantify the agreement of each 
quintessence cosmology with observational data, the cos- 
mological model that is most commonly referred to as the 
standard, concordance cosmological model with w = — 1 
erg = 0.9 (solid line, right panel) seems to be in conflict with 
the extant observational data. For w — — 1, lowering the 
normalization of the power spectrum to the value erg ~ 0.74 
(solid line, left panel) can greatly alleviate this discrepancy 
(Alam, Bullock, & Weinberg 2002; Zentner & Bullock 2002; 
McGaugh et al. 2003). However, a w = —0.5 model with the 
same as ~ 0.74 normalization (dotted line, left panel) does 
not do as well because earlier structure formation produces 
higher galactic central densities. 

Notice that the trends with A^ 2 for fixed V ma x do not 
scale as might be expected from the c v ir trends at fixed mass 
shown in Figure Hfjl This is because Ay/ 2 is a physical mea- 
sure of density and it increases not only with c v ir but also 
with A v i r (i.e. haloes are defined with respect to different 
overdensities). This effect is most apparent when comparing 
the right panels of Figures [TU1 and llll Though the concen- 
trations of haloes with w = —0.5 and as = 0.6 (dotted-line, 
right panel, Figure HOt are much lower than those in the 
standard w = —1.0, as = 0.9 case (solid line, right panel), 
the actual densities of those haloes are roughly the same in 
the right panel of Figure [TT1 This is because w = —0.5 mod- 
els have higher virial densities (see Figure 4). Even though 
haloes in the low-normalization w = —0.5 model tend to 
have the same (rapid-collapse) formation epoch as those in 
the higher- normalization w = — 1 model, the higher virial 
densities in the former model make the haloes denser over- 
all. 

By inspecting Figure llll we can immediately determine 
that models with w < — 1 and moderately low erg (erg ~ 0.7— 



0.8) can bring theoretical predictions to rough agreement 
with rotation curve data from low surface brightness and 
dwarf galaxies. Though it is clear that sufficiently decreasing 
erg can bring any model into accord with the median of the 
data, a w = —0.5 model would require erg < 0.6. Such a low 
normalization would be nearly impossible to reconcile with 
z — cluster abundance data. Thus from the standpoint of 
quintessence, models with w < — 1 seem mildly favored by 
galaxy density data. Conversely, models with w as high as 
~ —0.5 are strongly disfavored by galactic rotation curves 
coupled with only a weak prior on the normalization of the 
power spectrum. 

Note that none of the models can easily account for the 
very low data points. Nevertheless, the scatter in the data is 
not extremely large compared to the scatter expected from 
the halo-to-halo variations observed in N-body simulations. 
For example, at Vmax = 80 km/s, the lcr scatter in N-body 
simulations is cr(log(Av/ 2 )) — 0.37 while the lcr scatter in 
all 67 data points is <r(log(Avy2)) — 0.41. This suggests 
that lowering the median of the theoretically predicted cen- 
tral densities, perhaps by a reduction in erg or invoking a 
tilted or running power spectrum that reduces power on 
galaxy scales, or as we discuss here, by invoking w < — 1 
quintessence, may be sufficient to bring the predictions into 
good agreement with the data. Yet, we must bear in mind 
that our calculations are approximate. The most obvious 
omission is that all of our calculations are based on N-body 
simulations that contain no baryons. The effects of baryonic 
contraction are likely to be small in LSB galaxies (e.g., de 
Blok & McGaugh 1997) and would tend to drive rotation 
curves to higher values or Ay/ 2 and V max in the simplest 
models (Blumenthal et al. 1986). This serves only to increase 
the apparent discrepancy. Additionally, rotation curve mea- 
surements may yet be subject to poorly-understood system- 
atic effects in the reduction of the observational data (Swa- 
ters et al. 2003). Currently, it is difficult to draw a firm 
conclusion. 
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Figure 11. The central density parameter as a function of the maximum halo circular velocity. Dotted, solid, and dash-dot lines in each 
panel refer to w = —0.5,-1, and —1.5, respectively. The left panel is for a fixed <rg = 0.74 and the right panel is normalize to match the 
CMB: erg ~ 0.6 , 0.9, and 1.0 for w = —0.5,-1, and —1.5, respectively. The error bar in the right panel shows the expected theoretical 
1 — c scatter in Ay/2 due to the scatter in c v ; r . The points are for observed LSB and dwarf galaxies (see text for references). 



7 SUMMARY AND CONCLUSIONS 

Although the nature of dark energy is unknown, its effects on 
structure formation can be studied using numerical N-body 
simulations. We have performed a series of these simulations 
for a range of dark energy equation of state parameters. 
Confirming previous findings by Linder & Jenkins (2003), 
Klypin et al. (2003), Maccio et al. (2003), and Lokas et al. 
(2003) we show that the J01 formula provides a good fit to 
halo mass functions even in the presence of non-ACDM dark 
energy. We show that this is true for models with w < — 1 
as well. 

The density structure of dark matter haloes is also af- 
fected by dark energy. We have shown how the predictions 
of the B01 model are modified when dark energy with con- 
stant w is accounted for. As structure tends to form earlier 
in models with less negative w, halo concentrations tend to 
be somewhat higher in these models. These findings are in 
agreement with the results of Klypin et al. (2003) and qual- 
itatively agree with Dolag et al. (2003), although we probe 
a different range of masses than the latter. The larger num- 
ber of haloes with NFW profile fits and concentrations in our 
study allows us to quantitatively test the B01 model. We find 
that the original (F — 0.01, K = 4.0) over-predicts the con- 
centrations of haloes in our simulations by about ~ 12 — 15%. 
However, the shape of the mass-concentration relation that 
we find is the same as in B01, and we find that a slightly 
modified set of the B01 parameters (F = 0.01, K — 3.5) 
matches our haloes well. This offset may likely be caused 
by the lower force resolution of our GADGET simulations 
compared to the adaptive- refinement code ART used in BOT 
For the haloes in our simulations the adopted B01 model 
accurately reproduces the median concentration-mass re- 
lation over a range of masses from M v ir ~ 6 x 10 11 to 
M v i r ~4x 10 13 /2 -1 M Q . We confirm that for a fixed mass 
halo concentration decrease with redshift as 1/(1 + 2), at 
least out to z ~ 2.5. 

Interestingly, we find that halo concentrations are more 



easily understood when the halo virial radius is defined in 
terms of a cosmology-dependent virial overdensity rather 
than by one that uses a fixed overdensity of p/p C rit = 200. 
The result supports one of the (previously-untested) as- 
sumptions of the original B01 model. Specifically, we argue 
that halo densities in different cosmological models are in- 
fluenced both by changes in the overall virialization process 
of haloes as well as by changes in epoch when the halo cores 
collapse. As noted in it is important to include both 
of these physical processes when comparing predictions for 
galaxy densities to real data, as in Figure ITTI 

Having confirmed that the B01 model correctly de- 
scribes the scaling of halo concentrations as a function of 
mass and redshift even in cosmologies with w —1, we 
have investigated the effects of dark energy on a compari- 
son between model predictions and observations of central 
halo densities. Zentner and Bullock (2002) found that the 
observed distribution of A V /2 as a function of Vmax is in- 
consistent with the predictions of the B01 model for ACDM 
and erg = 0.9. The model predicts haloes that are simply 
too concentrated (see also Primack 2003). A lower value of 
<rg = 0.75, as used in our simulations, can alleviate this dis- 
crepancy, but such models may face other difficulties regard- 
ing early reionization (Somerville, Bullock, & Livio 2003) 
and possibly with reproducing the properties of halo sub- 
structure (Zentner & Bullock 2003). Including the effects 
of dark energy, we find that for models with w > — 1 the 
problem is exacerbated because haloes collapse earlier and 
because they have higher virial overdensities. Note that for 
the extreme case of w = —0.5, even a normalization as 
low as ag = 0.6 seems disfavored by the data. Thus one 
interesting conclusion is that the rotation curves of galax- 
ies coupled only with a weak prior on the normalization of 
the power spectrum of density fluctuations seem to disfa- 
vor quintessence models with w significantly larger than —1 
without measuring the expansion history of the Universe, as 
is done in SNIa analyses. Models with w < — 1 do better, 
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and can tolerate higher values of <xg ~ 0.8, but not high as 
high as as ~ 1, as is suggested the CMB normalization. 2 

As we have pointed out throughout this paper, there 
have been a number of previous studies of the effects of dark 
energy (m^ —1) on dark matter halo abundances and con- 
centrations. Here we summarize how this work distinguishes 
itself from these past efforts. We have shown for the first 
time how results for A v i r and the linear transfer function are 
modified for dark energy with w < — 1. In addition we have 
expanded on previous studies (Schuecker et al. 2003) of S c 
in phantom cosmologies. We have pointed out that there are 
significant differences in the amplitude and lu-dependence of 
as between normalizations based on the CMB and galaxy 
cluster abundances. For the first time we have been able 
to draw a distinction between a definition of concentration 
based on a cosmology-dependent R v i T and constant -R200, 
and found that our simulations prefer the former. Unlike 
previous studies of the effects of dark energy on halo con- 
centrations (Klypin et al. 2003; Dolag et al. 2003), which 
analysed a small number of pre-selected haloes, we have per- 
formed an analysis of ~ 1600 haloes in each simulation. This 
allowed us to directly test the B01 model for c(M, z) with 
statistical significance. Furthermore these previous studies 
focused on group- and cluster-sized haloes, whereas we have 
extended this study down to galaxy mass haloes. Lastly, we 
have shown for the first time how different values of w and 
power spectrum normalizations affect the theoretical predic- 
tions of Ay/ 2 vs. Vmax- A comparison to observations shows 
that models with w as low as -0.50 are strongly disfavored 
for any value of erg, standard ACDM models require as < 0.8, 
and models with w < — 1 can accomodate higher values of 
as- 

As future observations further constrain the nature of 
dark energy, it will become necessary to extend these types 
of studies to more realistic models of dark energy. Future 
simulations with higher mass and force resolution, includ- 
ing the effects of baryons, and incorporating more realistic 
dark energy models will further advance our understanding 
of the interplay between cosmology and dark matter halo 
structure. 
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2 Note that using K = 3.5 instead of K = 4.0 would not change 
any of the conclusions regarding galaxy rotation curves and the 
central density problem. 
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